# This script computes agency-level overexposure scores (without regressions)
# for Figure 3 of the main manuscript.

library(tidyverse)
library(lfe)
library(ggplot2)
load("combined_dta_long_483.rdata")

exp_weight = function(x, decay=0.00005, threshold = 483){
  y = exp(-decay*(x^2))
  y[x>threshold] = 0
  return(y)
}

dta = dta %>%
  mutate(
    w = exp_weight(d)
  ) %>%
  filter(year > 2009 & year < 2020)


all_dta_ori9 = data.frame(agency = NA, black_post = NA, black_arrest = NA, lrp=NA, type=NA)
crime_classes = unique(dta$crime_class)
crimes = unique(dta$crime)
q = length(crimes)
p = length(crime_classes)
for(i in c(1:q)){
  ct = crimes[i]
  cat(ct, "\n")
  ori9_post = dta %>%
    filter(crime==ct & type=='post') %>%
    group_by(agency) %>%
    summarize(
      black_post = sum((black>0)*w, na.rm=T) / sum(w)
    )
  
  ori9_arrest = dta %>%
    filter(crime==ct & type=='arrest') %>%
    group_by(agency) %>%
    summarize(
      black_arrest = sum((black>0)*w, na.rm=T) / sum(w)
    )
  
  new_results = full_join(ori9_post, ori9_arrest) %>%
    mutate(
      lrp = black_post - black_arrest
    )
  new_results$type = ct
  all_dta_ori9 = bind_rows(all_dta_ori9, new_results)
}

for(i in c(1:p)){
  ct = crime_classes[i]
  cat(ct, "\n")
  ori9_post = dta %>%
    filter(crime_class==ct & type=='post') %>%
    group_by(agency) %>%
    summarize(
      black_post = sum((black>0)*w, na.rm=T) / sum(w)
    )
  
  ori9_arrest = dta %>%
    filter(crime_class==ct & type=='arrest') %>%
    group_by(agency) %>%
    summarize(
      black_arrest = sum((black>0)*w, na.rm=T) / sum(w)
    )
  
  new_results = full_join(ori9_post, ori9_arrest) %>%
    mutate(
      lrp = black_post - black_arrest
    )
  new_results$type = ct
  all_dta_ori9 = bind_rows(all_dta_ori9, new_results)
}

cat("all", "\n")
ori9_post = dta %>%
  filter(type=='post') %>%
  group_by(agency) %>%
  summarize(
    black_post = sum((black>0)*w, na.rm=T) / sum(w)
  )

ori9_arrest = dta %>%
  filter(type=='arrest') %>%
  group_by(agency) %>%
  summarize(
    black_arrest = sum((black>0)*w, na.rm=T) / sum(w)
  )

new_results = full_join(ori9_post, ori9_arrest) %>%
  mutate(
    lrp = black_post - black_arrest
  )
new_results$type = "all"



all_dta_ori9 = bind_rows(all_dta_ori9, new_results)
all_dta_ori9 = all_dta_ori9[-1,]
write_csv(all_dta_ori9, "all_exposure_by_ori9.csv")